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Synopsis 

This  project  falls  ■within  the  realm  of  the  USAF  mission  of  developing  air  defense  systems, 
a  component  of  which  is  the  ongoing  Airborne  Laser  (ABL)  Program.  The  recent  congres¬ 
sional  mandate  to  develop  an  airborne  missile  defense  system  for  the  US  will  bolster  the  ABL 
program,  since  laser  technology  is  a  viable  tool  for  such  an  effort.  One  of  the  major  challenges 
of  the  ABL  program  is  the  development  of  laser-beam  propagation  codes.  Such  codes  must 
account  for  the  beam  propagation  through  an  extended  turbulent  medium  of  the  bulk  of 
the  troposphere  and  stratosphere.  From  the  operational  perspective  (Tactical  Decision  Aid 
Support),  rigorous  modeling  of  the  refractive  index  structure  function  for  long  horizontal  or 
nearly  horizontal  path  under  high  scintillation  conditions  characteristic  of  the  atmosphere 
is  imperative.  In  spite  of  demonstrated  importance  of  stratification  on  atmospheric  optical 
turbulence,  currently  available  ABL  phase  screen  theories  hinge  upon  isotropic  Kolmogorov 
spectrum  (Rytov’s  theory),  and  hence  cannot  represent  large  amplitude  fluctuations  of  at- 
mosphric  stratified  turbulence.  In  the  context  of  such  turbulence  the  collusion  between  the 
stratification  and  shear  leads  to  many  intriguing  phenomena  such  as  the  fiormation  of  thin, 
elongated  turbulent  layers  (pancakes)  and  instabilities  (such  as  Kelvin-Helmholtz  (K-H)  bil¬ 
lowing)  that  ultimately  break  down  into  turbulence.  The  turbulence  so  generated  is  often 
patchy  and  temporally  intermittent,  characterized  by  strong  anisotropy.  It  produces  strong 
optical  scintillation  due  to  refractive  index  fluctuations,  which  needs  to  to  be  quantified 
accurately  in  developing  beam  control  concepts  for  atmospheric  laser-beam  propagation. 

In  this  project,  direct  numerical  simulations  have  been  performed  to  study  the  dynamics 
of  an  inhomogeneous  stratified  shear  flow  that  models  an  atmospheric  jet  in  the  tropopause 
in  the  ABL  context.  The  basic  state  is  characterized  by  a  jet  centered  at  the  tropopause 
in  which  the  density  stratification  is  vertically  non-uniform.  Small  to  moderate  background 
stratifications  are  selected,  and  simulations  are  conducted  for  a  range  of  Reynolds  and  Froude 
numbers.  A  spectral  domain  decomposition  method  that  is  particularly  suitable  for  simu¬ 
lations  of  nonuniformly  stratified  stratified  shear  flows  is  developed  to  generate  the  desired 
turbulent  jet,  and  quasi-equilibrium  flow-fields  are  obtained  by  long-time  integration  of  gov¬ 
erning  equations.  The  structure  of  the  mean  flow  and  turbulence  fields  are  calculated,  which 
are  interpreted  using  relevant  length  scales  (Ozmidov,  buoyancy,  shear,  Ellison)  and  Richard¬ 
son  number  profiles.  The  ratios  of  the  Ellison  to  buoyancy  scales  are  much  smaller  than  unity 
at  the  jet  core  and  approach  unity  at  the  edges,  confirming  that  mechanical  turbulence  pre¬ 
vails  in  the  jet  core,  whereas  nonlinear  waves  and  stratification  effects  are  significant  at  the 
edges.  The  jet  core  is  found  to  support  sustained  mechanical  (active)  turbulence,  outside  of 
which  lay  a  region  of  patchy  turbulence  and  non-linear  gravity  wave  activity  characterized 
by  spatially  decaying  velocity  fluctuations  and  strong  temperature  fluctuations.  Detailed 
energy  budgets  show  how  energy  is  partitioned  within  the  flow,  including  the  transport  of 


1 


« 


energy  from  the  jet  to  its  immediate  vicinity  by  non-linear  gravity  waves.  Applications  in¬ 
clude  studies  of  the  vertical  variability  of  the  optical  refractive  index  function  Cl,  which  is 
of  primary  relevance  to  ABL  propagation  codes. 

The  following  personnel  has  been  supported  by  the  grant: 

•  Prof.  Alex  Mahalov,  Department  of  Mathematics  and  Mechanical  k  Aerospace  Engi¬ 
neering,  Center  for  Environmental  Fluid  Dynamics,  Arizona  State  University 

•  Prof.  Basil  Nicolaenko,  Department  of  Mathematics,  and  Center  for  Environmental 
Fluid  Dynamics,  Arizona  State  University 

•  Prof.  HJS  Fernando,  Center  for  Environmental  Fluid  Dynamics,  Arizona  State  Uni¬ 
versity 

•  Dr.  Binson  Joseph,  Post-Doctoral  Fellow,  Arizona  State  University 

•  Dr.  Frank  Tse,  Post-Doctoral  Fellow,  Arizona  State  University 

•  Dr.  Rafael  Pacheco,  Post-Doctoral  Fellow,  Arizona  State  University 

•  Prof.  A.  Babin,  Collaborator,  University  of  California,  Irvine 

•  Doctoral  students  Markus  Trahe  and  Bong-Sik  Kim 

Detailed  Technical  Report 
1  Introduction 

Some  of  the  challenges  in  Air  Borne  Laser  (ABL)  performance  predictions  are  tied  to  the  ef¬ 
fectiveness  of  propagation  codes.  The  latter  typically  use  discrete  grids  called  phase  screens 
to  atmospherically  induce  phases  to  be  simulated.  Key  questions  in  this  process  concern 
the  number  and  size  of  screens,  the  grid  spacing  and  the  details  of  using  atmospheric  data 
and/or  models  to  generate  the  field  strengths  at  the  grid  points.  Atmospheric  variability  gen¬ 
erates  spatial  intermittency,  non-stationarity  in  turbulence,  which  impacts  on  the  variability 
of  statistical  characteristics  of  C®,  the  refractive  index  structure  function.  A  fundamental 
characteristic  of  the  middle  atmosphere  (tropopause  and  stratosphere)  is  stratification  and 
layering:  optical  turbulence  occurs  in  thin  layers  (pancakes  or  “Mini”).  The  tropopause, 
by  definition,  is  a  region  of  sharp  vertical  temperature  gradients  at  and  above  which  strat¬ 
ification  and  layering  are  important;  its  altitudes  vary  from  around  12  km  at  midlatitudes 
to  8  km  in  polar  areas.  The  vertical  extent  of  pancakes  varies  from  100-150  m  to  30  m 
higher  in  the  stratosphere.  Balloon  data  show  persistent,  near  constant  altitude,  layers  of 
nearly  100  km  horizontal  extent.  The  extensive  variability  of  the  structure  of  the  tropopause 
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considerably  impacts  on  the  structure  of  turbulence.  There  is  a  demonsrated  need  to  better 
model  atmospheric  turbulence  along  long  quasi-horizontal  paths. 

In  the  atmosphere,  turbulence  is  generated  by  either  convective  or  shear  instabilities.  In 
the  middle  atmosphere  which  is  characterized  by  stable  background  stratification  and  shear, 
these  instabilities  and  turbulence  arise  primarily  from  two  sources: 

•  direct  instabilities  (such  as  shear  at  jet  streams,  convective  plumes); 

•  gravity  waves  induced  instabilities. 

In  the  middle  atmosphere  stable  background  stratification  and  shear  tend  to  influence 
motion  scales  that  are  larger  than  certain  :outef  scales  (e.g.  Ozmidov,  buoyancy,  Ellison). 
The  scales  of  active  Kolmogorov  turbulence  are  those  scales  which  are  greater  than  the  Kol¬ 
mogorov  scale  but  smaller  than  any  outer  scale.  The  Ozmidov  scale  is  regarded  as  the  outer 
scale  above  which  the  fluid  turbulence  is  significantly  affected  by  buoyancy.  Also,  the  dy¬ 
namics  of  a  stably  stratified  turbulent  flow  depend  on  the  rate  at  which  fluid  elements  mix 
each  other  and  change  their  density  in  the  interior  of  the  flow.  The  vertical  displacements  of 
the  fluid  elements  are  constrained  by  the  amount  of  energy  in  the  turbulence  to  lie  within  a 
buoyancy  scale  distance  of  the  order  of  the  buoyancy  scale  Lb.  If  the  vertical  velocity  fluc¬ 
tuations  due  to  internal  waves  are  large  compared  to  those  due  to  turbulence,  the  buoyancy 
scale  Lb  is  larger  than  L0  and  vice  versa.  The  scale  Le  known  as  the  overturning  or  Ellison 
scale  uses  the  density  fluctuations  as  an  indicator  of  turbulent  motions.  The  shear  length 
scale  Ls  represents  the  lower  bound  of  eddy  sizes  that  are  deformed  by  background  shear,  Lb 
is  the  buoyancy  scale  defined  above  and  the  scale  Ld  can  be  interpreted  as  the  integral  scale 
of  embedded  turbulence.  When  Ld  >  Ls  turbulence  interacts  with  and  extracts  energy  from 
the  background  shear,  thus  allowing  the  turbulent  intensity  to  grow.  The  stratification  can 
also  interact  with  turbulence  and  shear,  depending  on  the  relative  magnitudes  of  the  scales 
Ls,  Lb  and  Ld.  For  example,  when  Ld  <  Ls  <  Lb  the  shear  is  expected  to  interact  with 
stratification,  thus  producing  instabilities  (K-H  billowing),  gravity  waves  and  intermittent 
(patchy)  turbulence.  When  Lb  <  Ld  <  Ls  neither  the  instabilities  nor  shear-turbulence 
interactions  may  produce  turbulence,  and  the  existing  turbulence  in  the  flow  is  expected  to 
decay. 

Atmospheric  variability  generates  spatial  intermittency,  non-stationarity  in  turbulence, 
which  impacts  on  the  variability  of  statistical  characteristics  of  C2,  the  refractive  index 
structure  function.  For  classical  homogeneous  isotropic  Kolmogorov  turbulence,  it  is  defined 
by  Dn(r)  =  C^r2^,  where  Dn(r )  Is  the  refractive  index  structure  function  in  the  inertial 
convective  range: 

Dn(r)  =<  (n(x)  -  n(x  +  r))2  >, 

where  the  bracket  <,  >  denotes  ensemble  averaging.  The  three-dimensional,  nonstationary, 
intermittent,  anisotropic  structure  of  turbulence  in  the  atmospheric  field  is  coupled  with 
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optical  turbulence  in  thin  layers  (pancakes).  Balloon  measurements  have  been  used  to  probe 
into  the  vertical  structures  of  the  atmosphere.  Observations  so  obtained  demonstrate  the 
considerable  impact  of  sharp  vertical  temperature  gradients  in  the  tropopause  and  the  lower 
stratosphere  and  the  extensive  variability  of  the  structure  of  the  tropopause. 

The  following  are  some  major  issues  regarding  atmospheric  turbulence  impacting  on  ABL: 

•  Horizontal  pancake  turbulence  has  a  very  important  impact  on  C\\  very  little  is  exper¬ 
imentally  known  about  the  horizontal  structure  of  turbulence  and  C\\ 

•  Very  simple  models  such  as  CLEAR-1  and  “onion  skin”  perform  better  than  Rytov’s 
theory  models  but  are  still  not  adequate  for  turbulence  realizations  on  phase  screens  in 
propagation  codes; 

•  One  needs  3D  slew  model  that  includes  variability  due  to  intermittency,  nonstationarity 
and  stratification  (vary  C\  appropriately); 

•  Modeling  anisotropic  C\  along  paths  spanning  whole  pancakes  as  well  as  the  strong 
shearing  domains  between  pancakes  is  an  intricate  issue  that  need  to  be  investigated  by 
probing  into  the  structure  of  pancakes; 

•  Weather  front  dynamics  also  impact  via  intense  velocity  shearing  on  strong  variability 

of  Cl 

In  this  project  we  have  focused  on  investigating  turbulence  in  a  non-homogeneous  strat¬ 
ified  shear  flow  where  the  velocity  profile  takes  the  form  of  a  jet  (cf.  figure  1).  In  the  upper 
troposphere  and  stratosphere  localized  regions  of  three-dimensional  turbulence  arise  through 
shear  instability  or  through  breaking  of  inertia-gravity  waves.  Upper-level  atmospheric  jet 
regions  are  also  strongly  influenced  by  stable  background  stratification  (cf.  figure  2).  The 
flow  belongs  to  the  class  of  nonuniformly  stratified  (with  a  doubling  of  buoyancy  frequency 
across  the  jet)  shear  flows,  whose  turbulence  characteristics  are  not  well-understood.  We  set 
up  our  model  based  on  jet  streams  encountered  in  the  Earth  atmosphere  at  the  tropopause. 
This  region  corresponds  to  the  transition  between  the  tropopause  and  the  stratosphere  which 
one  observes  at  an  altitude  about  10  to  15  km  depending  on  the  latitude.  We  consider  a 
nonuniformly  stratified  jet  centered  at  the  vertical  coordinate  z  —  0.  The  flow  consists  of  a 
jet  core  surrounded  by  two  shear  layers;  the  layer  above  has  a  negative  shear  and  a  stronger 
temperature  gradient,  while  the  layer  below  has  a  positive  shear  and  presents  a  weaker  tem¬ 
perature  gradient,  with  the  buoyancy  (Brunt-Vaisala)  frequency  being  reduced  by  a  factor 
of  two.  This  latter  configuration  is  typical  of  the  jet  streams  in  the  tropopause  (cf.  figures 
1  and  2). 

The  gradient  Richardson  number,  which  quantifies  stratification  and  shear  effects,  is  low 
within  the  jet  core  (except  at  the  very  center)  and  velocity  fluctuations  are  maximal  there, 
thus  providing  sustained  turbulence  (cf.  figures  3a  and  3b).  It  increases  towards  the  jet  edges 
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where  the  effect  of  stratification  tends  to  reduce  turbulence  with  the  shear  length  scale  ex¬ 
ceeding  the  buoyancy  outer  scale.  Tropopausal  turbulence,  driven  by  shear  instabilities  on  ei¬ 
ther  side  of  the  jet  axis,  results  in  mean  stratification  with  a  notch  in  Brunt- Vaisala  frequency 
profile,  a  configuration  favoring  gravity  wave  emission.  The  flow  in  the  vicinity  of  the  edges 
of  the  jet  is  locally  out  of  equilibrium  and,  in  particular,  produces  nonlinear  gravity  waves 
which  travel  further  away  and  break  ([Sutherland  k  Peltier  (1995)]).  Regions  far  from  the  jet 
edges  have  weaker  velocity  fluctuations  than  in  the  core,  and  hence  meager  turbulent  mixing, 
although  potential  energy  and  temperature  fluctuations  remain  appreciable.  The  above  sce¬ 
narios  are  in  agreement  with  aircraft  observations  of  [Bedard,  Canavero  k  Einaudi  (1986)]. 
Critical  levels  where  waves  extract  energy  from  the  mean  flow  correspond  to  those  regions  of 
enhanced  turbulence,  where  the  Richardson  number  is  below  0.25,  hence  generating  Kelvin- 
Helmholtz  instabilities  ([Kaltenbach,  Gerz  k  Schumann  (1994)]). 

Our  focus  is  on  the  dynamics  and  vertical  variability  of  scales  that  are  poorly  resolved  in 
mesoscale  meteorological  codes  such  as  MM5  and  WRF  ([Grell,  Dudhia  k  Stauffer  (1995)]). 
Such  knowledge  is  important  not  only  for  parameterization  and  modeling  of  tropopausal  tur¬ 
bulence  and  associated  mixing  in  mesoscale  models  under  stably  stratified  conditions  (above 
the  atmospheric  boundary  layer),  but  also  for  determining  refractive  index  structure  func¬ 
tions  pertinent  to  electromagnetic  wave  propagation  across  the  tropopause  ([Beland  (1993)]; 
[Eaton  k  Nastrom  (1998)];  [Nastrom,  Gage  k  Ecklund  (1986)]; 

[Dalaudier,  Sidi,  Crochet  k  Vernin  (1994)]).  Vertical  scales  controlling  the  size  of ‘sheets’  in 
the  atmospheric  temperature  field  have  been  directly  evidenced  by  the  analysis  of  field  mea¬ 
surements;  the  main  dynamical  properties  of  such  layers,  whether  strongly  mixed  or  ‘calm’, 
can  be  characterized  not  only  by  the  gradient  Richardson  number  but  also  by  various  outer 
scales  of  turbulence  ([Alisse  k  Sidi  (2000)]).  We  investigate  how  the  dynamics  and  dominant 
physical  processes  in  the  stratified  jet  are  reflected  in  various  length  scales.  In  particular,  we 
study  vertical  variability  of  outer  scales  and  their  dependence  on  the  background  synoptic 
scale  jet  stream  and  the  background  Brunt-Vaisala  profile  across  the  tropopause.  Beyond 
a  sufficient  threshold  of  resolution  (especially  in  the  vertical  direction)  our  DNS  simula¬ 
tions  demonstrate  saturation  with  increased  numerical  resolution  of  the  vertical  variability 
curves  for  ratios  of  various  outer  scales;  cf.  figure  3b  for  Rig  =  L\jL\  and  figure  23b  for 
the  temperature  fluctuation  parameter  Le/Lb,  where  Ls  is  the  shear  scale,  Lh  is  the  buoy¬ 
ancy  scale  and  Le  is  the  Ellison  scale.  These  resolution  independent  ratios  are  an  effective 
gauge  of  the  variability  of  inhomogeneous  mixing  in  a  nonuniformly  stratified  tropopause  jet. 

Only  a  very  few  studies  have  been  reported  with  regard  to  three-dimensional,  high- 
resolution  numerical  simulations  that  employ  realistic  models  of  the  nonuniformly  stratified 
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tropopause  jet.  [Sutherland  k  Peltier  (1995)]  and  [Smyth  k  Mourn  (2002)]  implemented  a 
jet-like  velocity  profile  as  the  initial  condition;  the  Brunt-Vaisala  profiles  used  in  their  model 
was  also  doubled  from  the  lower  to  the  upper  domains.  However,  their  simulations  are  only 
two  dimensional  and  the  jet  profiles  tend  to  be  smeared  as  there  is  no  sustained  momen¬ 
tum  source  in  their  model,  which  caused  the  shear  to  steadily  decrease.  In  the  atmosphere, 
on  the  other  hand,  the  jet  is  maintained  quasi-steady  by  the  synoptic  scale  forcing,  which 
needs  to  be  considered  in  modeling.  It  should  be  noted  that  there  are  many  numerical 
studies  of  stratified  turbulent  shear  flows,  which  impose  a  constant  mean  shear  and  stratifi¬ 
cation  in  the  vertical  direction.  The  horizontal  directions  are  then  regarded  as  periodic  and 
the  vertical  direction  shear  periodic.  The  work  by  [Gerz,  Schumann  k  Elghobashi  (1989)]; 
[Holt,  Koseff  k  Ferziger  (1992)];  [Kaltenbach,  Gerz  k  Schumann  (1994)], 

[Carnevale,  Briscoline  &  Orlandi  (2001)]  and  [Jacobitz,  Sarkar  k  Van  Atta  (1997)]  belong  to 
this  category.  Turbulent  flows  with  mean  shear  and  stratification  are  typically  temporally 
evolving  and  hence  do  not  reach  stationary  states,  except  at  a  certain  ‘stationary’  Richard¬ 
son  number  for  which  the  production  is  balanced  by  the  dissipation.  Asymptotic  stationary 
states  can  be  attained,  however,  where  normalized  variables  (e.g.  ratios  of  individual  compo¬ 
nents  of  energy  to  total  energy)  reach  constant  values  after  several  turnover  periods.  Another 
approach  in  the  study  of  shear-stratified  turbulence  uses  decaying  simulations.  An  example 
is  the  DNS  studies  of  [Galmiche,  Thual  k  Bonneton  (2002)]  dealing  with  non-uniform  ver¬ 
tical  mean  shear  and  non-uniform  stratification,  where  mean  velocity  and  density  profiles 
are  allowed  to  evolve  with  time.  In  such  simulations,  however,  the  flow  does  not  reach  equi¬ 
librium  and  it  remains  unclear  whether  the  results  are  sensitive  to  the  generally  transient 
nature  of  the  flow. 

While  previous  numerical  studies  have  made  important  contribution  to  the  study  of  shear- 
stratified  turbulence,  they  do  not  reproduce  exact  middle  atmospheric  situation:  wherein  the 
turbulence  can  reach  some  quasi-equilibrium  state  and  where  the  turbulent  kinetic  energy 
(TKE)  budgets  are  nearly  balanced.  Furthermore,  if  homogeneous  boundary  conditions  are 
used  in  the  vertical  direction,  there  is  no  net  heat  and  momentum  transfer  out  of  the  box.  As 
we  demonstrate  in  Section  3,  the  fluxes  play  an  important  role  in  redistributing  energy  in  the 
vertical  direction.  Without  such  a  mechanism,  the  variability  of  the  atmospheric  turbulence 
around  the  jet  cannot  be  reproduced.  It  has  been  shown  in  [Sutherland  k  Peltier  (1995)] 
that  the  momentum  flux  created  from  the  jet  produces  internal  (i.e.  non-topographical) 
gravity  waves  which  In  turn  produce  a  drag  on  the  mean  velocity.  These  phenomena  can¬ 
not  be  reproduced  by  conventional  shear  periodic  numerical  experiments,  the  application  of 
which  to  atmospheric  tropopause  simulations  in  the  presence  of  a  jet  stream  is  thus  limited. 
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Our  model  of  the  atmospheric  tropopause  is  characterized  by  a  stratified  thin  layer  across 
which  the  buoyancy  frequency  jumps  approximately  by  a  factor  of  two  (figure  2).  In  our 
three-dimensional  direct  numerical  simulations  (DNS),  the  streamwise  velocity  profile  that 
characterizes  the  basic  state  of  the  model  jet  was  set  to  a  gaussian  profile,  thus  produc¬ 
ing  shear  at  different  vertical  levels  characteristic  of  a  jet  stream.  The  governing  equations 
are  solved  for  the  above  flow  configuration  using  DNS  based  on  a  spectral  domain  technique 
developed  specifically  for  flows  characterized  by  non-uniform  background  shear  and  stratifica¬ 
tion  ([Tse,  Mahalov,  Nicolaenko  &  Fernando  (2001)]).  Due  to  inhomogeneity  in  the  vertical 
direction,  periodic  boundary  conditions  are  not  used  therein;  previous  numerical  studies  have 
assumed  such  periodicity.  In  our  simulations,  non-linear  shear  and  stratification  profiles  ad¬ 
just  to  roughly  stationary  values,  around  which  the  potential  and  kinetic  energies  fluctuate. 
Quasi-equilibrium  solutions  are  obtained  following  long-time  integration  of  governing  equa¬ 
tions.  The  robustness  of  our  results  has  been  verified  by  DNS  at  doubled  resolution  (1024 
vertical  levels).  Thus,  physical  results  reported  here  and,  in  particular,  those  involving  ra¬ 
tios  of  outer  scales  are  resolution  independent.  We  investigate  the  effects  of  inhomogeneity 
of  shear  and  stratification,  especially  those  properties  that  are  distinct  from  homogeneous 
stratified  shear  flows. 

The  next  section  describes  the  computational  framework  for  non-homogeneous  stratified 
shear  flows.  Section  3  presents  the  results  of  numerical  simulations  and  discussions  of  various 
mean  quantities  and  variances  at  quasi-equilibrium.  Vertical  variability  of  natural  length 
scales  (Ozmidov,  buoyancy,  shear  and  Ellison)  and  Richardson  number  profiles  are  analyzed 
in  Section  4.  Turbulent  budget  equations  are  investigated  in  Section  5.  Conclusions  and 
issues  for  future  work  are  discussed  in  Section  6. 

2  Computational  framework  for  non-homogeneous  strat¬ 
ified  shear  flows 

2A  Description  of  numerical  experiments 

The  governing  equations  are  the  three-dimensional  incompressible  Navier-Stokes  equations 
for  the  velocity  U  and  temperature  0  under  the  Boussinesq  approximation: 


§  —  ■ 

dP  _2rr 

-_  +  „vC+nt,i 

(2.1) 

dV 

dP  , 

aF  +  u'w  = 

+  !/V2V, 
ay 

(2.2) 

dW  „ 

dP  2 

-m+vvw  = 

-^  +  uv2w  +  gp(e-eR), 

(2.3) 
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-^  +  u-ve  =  kv2q  +  ne ,  (2.4) 

v-u  =  0,  (2.5) 

where  U  =  (U,  V,  W)  are  the  three  components  of  the  velocity  in  the  streamwise,  spanwise 
and  vertical  directions  (denoted  as  x,y,z)  respectively,  QR  is  the  constant  reference  temper¬ 
ature;  P  is  the  pressure;  u,  k,  (3  and  g  are  the  molecular  viscosity,  molecular  diffusivity, 
thermal  expansion  coefficient  and  gravitational  acceleration,  respectively.  The  horizontally 
homogeneous  terms  II ^  and  11©  are  the  vertically  dependent  momentum  and  thermal  sources. 
The  two  source  terms  enable  to  obtain  a  basic  (unstable)  jet  stream  profile  and  a  desired 
vertically  variable  Brunt- Vaisala  frequency  profile.  They  represent  large  scale  momentum 
and  thermal  forcing  responsible  for  maintaining  the  jet. 


From  now  on,  the  total  instantaneous  variables  are  represented  by  upper  case  letters 
(U,  V ,  W,  0).  They  are  decomposed  into  two  parts  :  basic  state  and  perturbations.  The 
variables  in  the  basic  state  are  represented  by  upper  case  letters  with  subscripts  (Ub  and 
©h  +  ©h),  and  perturbations  by  lower  case  letters  (u,v,w,B).  The  perturbation  is  further 
decomposed  into  two  components:  a  mean  (horizontally  averaged)  part  denoted  by  (•)  and  a 
fluctuating  component  denoted  by  primes.  Thus,  the  variables  U,  V,  W,  ©  are  decomposed 
as  : 


U  =  Ub(z)+u(x,y,z,t)  =  Ub(z)  +  (u)(z,t)  +  u'(x,y,z,t)t  (2.6) 

V  =  v(x,y,z,t)  =  (v)(z,t)  +v'(x,y,z,t) ,  (2.7) 

W  =  w(x,  y,  z ,  t)  =  w'(x,  y,  z,  t) ,  (2.8) 

0  =  Qii  +  9b(z) +  9(x,y,z,t)  =  0R  + Qb(z)  +  (9)(z,t)  +  9'(x,y,z,t)  .  (2.9) 


We  have  (u')=(v,)=(vji)  =  0.  The  mean  (horizontally  averaged)  values  of  U  and  0  contain 
a  part  from  the  basic  state  and  a  part  from  the  perturbation  (e.g.  U  =  Ub  +  (u)),  while  the 
mean  of  V  satisfies  (V  =  (v)).  The  pressure  can  also  be  decomposed  as  : 

P  =  Pb{z)+p(x,y,z,t)  =  Pb(z)  +  (p){z,t)  +p'(x,y,z,t) .  (2.10) 


For  the  basic  state,  where  the  perturbations  u,  v,  w  and  9  are  zero  and  the  fields  depend 
only  on  z,  the  governing  equations  are  reduced  to  : 


n„  +  4^  =  o, 


dPh 


dz 2 

n©  =  — k 


d2Qb 


dz  “O  '“dz2' 

Here  Yiu  is  obtained  by  assigning  to  Ub  the  following  gaussian  form: 


(2.11) 

(2.12) 

(2.13) 
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where  <*1  is  the  stiffness  parameter  and  Lz  is  the  half  depth  of  the  computational  box.  The 
velocity  profile  Ub(z)  is  equal  to  £4(0)  at  the  center.  The  momentum  source  11^  and  the 
thermal  source  II©  are  obtained  from  eq.  (2.11)  and  from  eq.  (2.12),  respectively.  The  mo¬ 
mentum  source  is  intended  to  maintain  a  sustained  synoptic/planetary  jet  stream  in  the 
model  (figure  1).  The  thermal  source  maintains  the  mesoscale  doubling  of  the  buoyancy 
frequency  across  the  tropopause.  Thus,  we  are  resolving  a  microscale  box  centered  on  a 
synoptic-scale  jet  stream  at  the  tropopause,  for  times  much  shorter  than  time  scales  of  jet 
dynamics.  The  jet  profile  for  Ub  is  shown  in  figure  lb,  with  a  =  ax  =  16.  This  profile  models 
the  measured  wind  profile  of  an  atmospheric  jet  shown  in  figure  la. 

In  eq.  (2.12),  the  buoyancy  force  due  to  background  temperature  gradient  is  balanced  by 
the  basic  state  pressure,  the  thermal  source  is  balanced  by  the  diffusion  of  the  background 
temperature.  Substituting  eq.  (2.9),  (2.10)  and  (2.12)  into  eq.  (2.3)  -  (2.4),  the  equations 
for  vertical  velocity  and  temperature  can  then  be  expressed  as  : 

^  +  U-W  =  -^  +  g/39  +  vV2W,  (2.14) 

l^  +  U-W  =  -W^  +  kV20;  (2.15) 

where  /?  is  the  thermal  expansion  coefficient.  Next,  let  N2K2(z)  =  g/3dQb(z)/dz  and 
1?  =  g/39/N  where  N  is  a  constant  (iV-1  has  the  unit  of  time)  and  K(z)  is  a  verti¬ 
cal  profile  factor,  we  obtain  the  rescaled  vertical  velocity  and  temperature  equations  (e.g. 
[Herring  &  Metais  (1989)]): 

dW  dv 

— +  U-VW  =  4-  N$  +  uV2W ,  (2.16) 

dd 

— +  U-Vi?  =  -WNK2  +  kV2$  .  (2.17) 

The  variable  1?  has  the  unit  of  length/time  (not  to  be  confused  with  9  which  has  the  unit 

of  temperature).  In  the  above  equation  N  and  K(z)  are  related  to  the  usual  Brunt- Vaisala 
frequency  profile  N2(z)  —  ©ff  as: 

N2K\z)  =  gfl- ^  =  p ej^(z) .  (2.18) 

In  our  simulations,  the  non-dimensional  K(z)  profile  has  the  following  form: 

=  1  +  e-4  +  l  ■  (2-19) 

The  profile  K  (z)  corresponds  to  a  temperature  profile  which  increases  monotonieally  from 
the  bottom  to  the  tropopause  level,  experiences  a  stiff  transition  at  the  tropopause  which 
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depends  on  the  stiffness  parameter  a2,  then  further  increases  reaching  at  the  top  twice  the 
value  at  the  bottom.  This  dimensionless  profile  models  the  doubling  of  Brunt- Vaisala  fre¬ 
quency  at  the  tropopause,  as  observed  in  balloon  field  measurement  by  [Beland  (1993)], 
figure  2  (doubling  roughly  from  0.01s-1  to  0.02s-1).  Equations  (2.13)  and  (2.19)  introduce 
two  parameters  ax  and  a2,  which  specify  the  inverse  of  the  external  length  scales  of  the  sys¬ 
tem.  The  first  number  ax  controls  the  amount  of  shear  along  the  jet  edges  while  the  second 
parameter  a2  controls  the  stiffness  in  Brunt-Vaisala  profile.  In  our  simulations,  the  same 
values  are  used  for  ax  and  a2  (a  =  ax  =  a2  =  16).  The  edges  of  the  jet  are  then  located  in 
the  vicinity  of  the  normalized  vertical  levels  zam±l. 


The  variables  can  be  non-dimensionalized  by  choosing  a  suitable  velocity  scale  UQ  =  £4(0) 
and  length  scale  La.  We  define  La  =  1  /a  and  z  =  az  =  z/La.  Hence,  in  terms  of  the  rescaled 
z  and  with  a  =  16,  the  computational  box  corresponds  to  -5.0  x  16  <  z  <  5.0  x  16.  Similarly, 
we  rescale  x  =  ax,  y  =  ay,  with  0  <  x ,  y  <  8x.  The  nondimensionalized  equations  then 
become  : 


dU  ~ 

af+u'w  = 

(2.20) 

o 

< 

II 

+  1  V2v , 
dy^Reo  ’ 

(2.21) 

dW  ~ 

af_+u'w  = 

+  -^v2w  +  , 

az  Re0  Fr 

(2.22) 

||  +  U-Vtf  = 
at 

—^WK2*^1  V2tf. 

Fr  Re$Pr 

(2.23) 

leaving  only  three  dimensionless  parameter  groups  in  the  equations:  the  Froude  number 
Ft  =  Ua/NLa,  initial  Reynolds  number  Re0  =  UaLa/v  and  Prandtl  number  Pr  =  v/k. 
Similarly,  the  source  now  has  the  form  : 


1  32£4 
Re0  dz2  ’ 


(2.24) 


The  governing  equations  (2.1),  (2.2),  (2.16)  and  (2.17)  are  solved  using  numerical  meth¬ 
ods  described  in  Section  2.2.  The  parameters  for  this  problem  are  k,  v,  £4(0),  N,  and  a; 
their  values  used  in  our  simulations  are  listed  in  table  1.  Case  1  is  regarded  as  the  reference 
case.  The  viscous  and  thermal  diffusivities  are  set  to  1  x  10-4  and  1.4  x  10-4,  respectively, 
for  cases  1-5  and  6  x  10-5  and  8.6  x  10-5  for  case  6,  resulting  in  a  Prandtl  number  v/k  of 
0.7  in  all  cases  as  for  air.  Case  2  has  higher  and  case  3  has  lower  stratification  compare  to 


10 


Table  1:  Physical  parameters  used  in  simulations  with  v  —  1  x  10~4,  k  =  1.4x  10-4  (cases 
1-5)  and  y  =  6x  1CT5,  «  =  8.6  x  10-5  (case  6).  Pr  —  0.7  in  all  cases. 


Case 

£4(0) 

resolution 

Re°  =  H 

N 

1 

4 

1282  x  512 

2500 

0.2 

2 

4 

1282  x  512 

2500 

1.0 

3 

4 

1282  x  512 

2500 

0.05 

4 

10 

2562  x  512 

6250 

0.2 

5 

20 

2562  x  512 

12500 

0.2 

6 

20 

5122  x  1024 

20833 

0.2 

case  1.  The  dependence  on  stratification  can  be  studied  through  the  first  three  cases.  Cases 
4,  5  and  6  have  higher  forcing  respectively,  as  evident  of  their  respective  Reynolds  numbers. 
Cases  4  and  5  are  simulated  with  a  resolution  of  2562  x  512  while  cases  1-3  have  a  resolution 
of  1282  x  512.  The  largest  Reynolds  number  simulation  in  case  6  has  the  highest  resolution 
of  5122  x  1024. 

The  vertical  gradient  of  Ub(z)  is  large  in  the  middle  region  (except  at  z  =  0)  and  corre¬ 
spondingly  Rig  =  (. NK(z))2/(dUb/dzf  is  small  there.  The  base  jet  profile  is  linearly  unstable 
(Rig  <  0.25).  The  shear  is  sufficiently  large  to  overcome  the  stability  due  to  stable  stratifi¬ 
cation  and  the  flow  develops  instabilities  and  turbulence.  Based  on  the  centerline  velocity 
of  the  jet  at  the  basic  state  £4(0)  and  a-1,  the  initial  Reynolds  number  is  given  in  table 
1.  Case  1,  2  &  3  have  the  same  Reynolds  number  but  different  background  stratification. 
Increasing  the  forcing  parameter  has  the  effect  of  increasing  the  Reynolds  number  (cases  4, 
5  and  6).  For  the  quasi-equilibrium  turbulent  state  we  define  turbulent  length  and  velocity 
scales  as  follows.  The  lengthscale,  Ld  =  q3/e  is  based  on  the  turbulent  kinetic  energy  (TKE) 
dissipation  e  and  q  ([Batchelor  (1953)]): 

e  =  uisijSij ) ,  q2  =  (u2  +  v2  +  w 2) ,  (2.26) 

where  Sy  =  1/2 (dui/dxj  +  duj/dxi )  is  the  rate-of-strain  tensor  ([Pope  (2000)])  and  ut  and 
Xi  (i  =  1,2,3)  correspond  to  our  u,  v,  w  and  x,  y,  z  notation  respectively.  The  values  of 
turbulent  Reynolds  number  Red  and  Froude  number  Frd  with  this  set  of  scales  are  given 
in  table  2  for  the  vertical  level  az  =  1.0.  The  scales  and  the  Reynolds  number  Red  vary 
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Table  2:  Values  of  parameter  groups  at  quasi-equilibrium  at  the  vertical  level  az  =  1.0; 
Froude  Frd  —  and  Reynolds  Red  =  ^  numbers  are  based  on  the  length  scale  Ld. 


Case 

1 

2 

3 

4 

5 

6 

Frd 

1.26 

0.41 

5.21 

1.98 

2.54 

1.11 

Red 

826 

534 

798 

1392 

2392 

4016 

with  vertical  levels,  being  highest  near  the  edges  of  the  jet  and  decay  rapidly  away  from  the 
edges;  see  figure  4  for  case  6  with  512  x  512  x  1024  resolution. 


2B  Numerical  methods 

Most  of  the  previous  numerical  simulations  of  shear-stratified  turbulence  have  been  in  the 
context  of  uniform  background  stratification  or  uniform  shear.  In  the  present  study  we 
use  a  spectral  domain  decomposition  method,  which  is  particularly  suitable  for  simulation 
of  flows  with  non-uniform  background  stratification  and  shear.  The  flow  is  assumed  to 
be  homogeneous  only  in  the  horizontal  directions  where  periodic  boundary  conditions  are 
used.  For  each  horizontal  wavenumber,  the  vertical  domain  is  then  broken  down  into  several 
subdomains  (figure  5).  Each  subdomain  is  separately  mapped  to  a  domain  {—\<zf<l}. 
For  non-overlapping  subdomains,  there  are  basically  two  solution  methods:  collocation  and 
variational  methods.  In  the  collocation  method,  the  variable  (e.g.  A)  within  each  subdomain 
j  is  interpolated  as  : 

XJ(Z')  =  f24Mz')i  -1  <  z'  <  1 ,  (2.27) 

i- 0 

where  n  is  the  order  of  interpolation,  and  the  functions  0i(z')  are  the  Lagrange  interpolants. 
The  local  coordinate  z'  is  chosen  to  be  the  Gauss  Lobatto  Legendre  points. 

The  differentiation  of  variables  are  carried  out  by  differentiating  the  interpolants.  Effi¬ 
cient  routine  for  the  differentiation  exists,  e.g.  [Fornberg  (1996)].  Domain  decomposition 
method,  however,  requires  patching  of  subdomain  boundaries  when  forming  the  global  dif¬ 
ferentiation  matrix.  In  inverting  the  2nd  order  differentiation  matrix  (e.g,  in  calculating  the 
pressure  Poisson  equation),  continuity  of  the  1st  order  derivative  is  required  as  subdomain 
boundary  condition.  This  is  usually  done  by  replacing  n  —  1  rows  (corresponding  to  n  —  1 
subdomain  boundaries)  of  the  global  matrix  by  the  1st  order  derivative  equations.  This 
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method,  however,  may  lead  to  numerical  disturbances  at  the  subdomain  boundaries,  partic¬ 
ularly  if  the  resolution  of  the  problem  is  marginal.  To  overcome  this  problem,  an  auxiliary 
set  of  subdomains  is  defined,  with  each  subdomain  in  the  auxiliary  set  strides  across  the  sub- 
domains  in  the  primary  set.  The  arrangement  is  shown  schematically  in  figure  5.  Another 
global  differentiation  matrix  is  formed.  Differentiation  or  integration  are  then  carried  out 
by  multiplication  or  inversion  of  a  combination  of  these  two  matrices.  When  this  method 
is  used  in  advection-diffusion  equation,  exponential  convergence  of  accuracy  is  observed  for 
increasing  order  of  interpolation.  The  method  has  been  tested  for  the  turbulent  channel 
flow  problems,  and  excellent  performance  has  been  noted.  Further  details  on  the  numerical 
method  are  given  in  [Tse,  Mahalov,  Nicolaenko  &  Fernando  (2001)]. 

The  time  discretization  follows  the  usual  pressure  projection  method,  which  separates 
the  time  derivatives  into  three  substeps.  The  nonlinear  terms  and  buoyancy  terms  are 
advanced  in  the  first  substep;  pressure  Poisson  equations  are  then  solved  directly  and  the 
pressure  terms  are  advanced  in  the  second  substep  while  the  dissipation  term  is  considered 
in  the  third.  Second  order  Adams-Bashforth  scheme  is  used  for  the  nonlinear  terms.  The 
calculation  of  derivatives  involved  in  the  nonlinear  terms  is  carried  out  in  physical  space  using 
collocation  methods.  The  third  substep  requires  calculating  the  viscous  terms  implicitly.  A 
total  of  eight  forward  FFT  transforms  and  four  backward  FFT  transforms  is  required.  In 
the  vertical  direction,  absorption  layers  are  added  to  the  top  and  bottom  boundaries  to 
remove,  or  at  least  reduce,  reflection  of  waves.  The  length  of  the  computational  box  shown 
in  figure  5b  is  7t/2  in  horizontal  directions  and  from  -5.0  to  +5.0  in  the  vertical  direction. 
Grid  resolution  in  our  work  includes  128  x  128  x  512,  256  x  256  x  512  and  512  x  512  x  1024 
(table  1).  There  are  127  (for  cases  1-5)  and  255  (for  case  6)  subdomains  vertically.  The 
interpolation  for  the  central  three  subdomains  is  5th  order  while  for  the  rest  is  4th  order. 
The  total  number  of  grid  points  in  the  vertical  is  then  power  of  two,  avoiding  one  processor 
carrying  an  extra  plane  for  computation.  The  width  of  the  subdomains  between  -0.062  to 
0.062  is  0.004,  and  it  gradually  increases  to  0.2  at  the  outer  region.  The  jet  stream  is  placed 
in  the  center  of  the  computational  domain.  This  region  requires  more  resolution  due  to 
high  shear  (figure  5).  The  program  is  parallelized  by  the  transposition  method  using  MPI. 
Typically,  for  resolution  of  128  x  128  x  512,  it  takes  15  seconds  to  advance  one  time  step 
on  a  cluster  of  32  SGI  R10000  processors  and  for  512  x  512  x  1024  resolution,  it  takes  »  5 
minutes  on  128  processors.  The  computations  are  done  on  the  Nirvana  massively  parallel 
supercomputer  at  Los  Alamos  National  Laboratory  and  the  DoD  ARL  MSRC. 
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2C  Quasi-equilibrium  state 


Long  time  integration  of  the  governing  eqs.  (2.1)  -  (2.5)  is  affected  by  the  presence  of  adia¬ 
batic  invariants  for  the  long  time  dynamics.  In  particular,  it  has  been  shown  in 
[Babin,  Mahalov  k  Nicolaenko  (1998)]  that  in  the  asymptotic  state  (after  several  periods  of 
oscillations  associated  with  wave  motions;  [Mahalov,  Nicolaenko  k  Zhou  (1998)])  the  hor¬ 
izontally  averaged  temperature  (d)  is  a  near  adiabatic  invariant  with  onl  y  a  small  drift 
in  time;  this  drift  is  confirmed  by  figure  6  in  our  DNS.  We  recall  the  mean  (horizontally 
averaged)  temperature  equations: 


d(d)  d^u/f)  _  d2{d) 
dt  dz  K  dz2 


(2.28) 


A  state  of  absolute  equilibrium  is  one  where  the  mean  temperature  (tf)  as  well  as  other 
(horizontally)  averaged  quantities  are  no  longer  changing  with  time.  We  define  quasi¬ 
equilibrium  to  be  a  state  where  (t?)  satisfies  eq.  (2.28)  with  a  small  but  finite  time  derivative 
(slow  steady  drift  as  in  figure  6);  and  where  all  other  means  are  quasi-stationary  and  the 
integrated  statistics  of  variances  show  little  fluctuation  (<  5%).  This  is  illustrated  in  figure 
7  for  a  typical  case  where  the  vertically  integrated  values  of  velocity  variances,  J(u'2)dz  (I), 
f(v'2)dz  (II)  and  f(w'2)dz  (III),  temperature  variance  f(d'2)dz  (IV)  and  the  rate  of  change  of 
horizontal  velocity  variance,  /  \d{u2)/dt\dz  (V),  are  plotted  against  the  large  eddy  overturn¬ 
ing  time  (S6-1  =  (max  \dUb/dz\)~l  ss  0.018),  Our  range  for  Sbt  is  longer  than  those  reported 
in  homogeneous  simulations  since  we  are  using  the  maximum  value  of  shear  at  basic  state. 


The  actual  shear  increases  from  zero  to  the  quasi-equilibrium  values  and  also  depends  on 
the  vertical  level.  All  the  curves  initially  overshoot  and  then  settle  down  to  a  relatively 
constant  value.  The  first  stage  of  evolution  corresponds  to  the  onset  of  Kelvin-Helmhotz 


instabilities.  Secondary  instabilities  then  kick-in,  with  transition  to  3D  turbulence  and  the 


values  of  the  variances  drop  and  begin  to  saturate.  The  same  trend  has  been  observed,  for 
example,  in  [Smyth  &  Mourn  (2000)].  After  the  velocity  variances  become  quasi-stationary, 
the  horizontally  averaged  mean  temperature  is  still  evolving,  albeit  slowly,  as  shown  in  figure 
6  for  case  1.  Physically,  quasi-equilibrium  states  can  be  obtained  if  the  synoptic  scale  driv¬ 
ing  force  has  a  much  longer  time  scale  (e.g.  see  figure  3  of  [Cullen  (2002)]).  The  turbulence 
so  obtained  is  in  local  equilibrium  with  TKE  production  in  large  scales  which  dissipates 
at  smaller  scales.  This  is  evident  from  the  spectra  in  quasi-equilibrium  state  presented  in 
[Tse,  Mahalov,  Nicolaenko  k  Fernando  (2001)],  which  span  several  decades. 


3  The  mean  profiles,  variances  and  covariances 

In  this  section,  based  on  our  long-time  DNS,  we  present  vertical  profiles  at  quasi-equilibrium 
for  the  mean  temperature,  three  velocity  components  and  the  corresponding  variances  and  co- 


14 


variances.  The  mean  (horizontal  averaged)  velocity  in  the  streamwise  direction  U  =  Ub  +  {u) 
at  quasi-equilibrium  is  shown  in  figure  8  for  the  cases  from  table  1.  Those  figures  are 
time-averaged  for  one  large  eddy  turnover  time  defined  by  2tt /(aq).  The  vertical  axis  is 
nondimensionalized  by  1/a  and  the  velocity  is  non-dimensionalized  by  17(0).  Hence  the  nor¬ 
malized  values  are  all  equal  to  unity  at  the  center.  Cases  1,  2  and  3  are  shown  in  figure  8a 
while  cases  4,  5  and  6  are  shown  in  figure  8b.  The  profiles  for  cases  5  and  6  are  essentially 
identical.  The  amount  of  stratification  does  not  have  much  effect  on  the  width  of  the  mean 
flow,  which  is  dictated  by  the  profile  of  the  momentum  source  11^.  There  is  a  counter  flow 
away  from  the  central  region.  Note  that  the  asymmetry  generated  by  the  doubling  of  the 
Brunt-Vaisala  profile,  is  more  prominent  in  case  2  than  in  the  other  cases.  The  peak  value 
of  U  =  Ub  +  (u)  in  this  case  is  also  larger  than  in  other  cases.  In  the  decaying  simulation 
by  [Galmiche,  Thual  k  Bonneton  (2002)],  an  increase  of  stratification  has  been  observed  to 
cause  transfer  of  energy  from  turbulence  to  the  mean.  Apparently,  we  encounter  a  similar 
situation  in  our  quasi-equilibrium  simulations. 

The  mean  temperature  profiles  ($),  which  is  the  mean  of  the  perturbation  of  the  basic 
profile,  is  a  slowly  evolving  quantity  that  can  serve  as  an  indicator  of  a  particular  quasi¬ 
equilibrium  state  at  that  vertical  level.  Their  graphs  are  plotted  in  figure  9,  with  values 
normalized  by  17(0)  (t?  has  the  unit  of  velocity).  The  values  are  negative  in  the  upper  do¬ 
main  and  positive  in  the  lower  domain.  As  will  be  seen  later  in  figure  20,  the  vertical  heat 
flux  is  negative  in  the  central  region  and  it  is  very  small  elsewhere.  As  a  result,  the  tem¬ 
perature  decreases  at  the  top  and  increases  at  the  bottom  due  to  buoyancy  flux  crossing  the 
center.  In  figure  9  the  peaks  for  case  2  are  the  largest,  followed  by  case  1  while  case  3  has 
the  smallest  peaks.  One  prominent  feature  is  that,  as  the  stratification  decreases,  the  peaks 
move  away  from  the  core  of  the  flow.  Consider  case  3  for  example,  (i?)  peaks  at  za  «  5.36 
and  -4.96,  at  distances  significantly  away  from  the  sides  of  the  jet  at  za  «  ±1.  There  are 
no  significant  differences  in  the  location  of  peaks  for  cases  4,  5  and  6. 

Figure  10a  shows  profiles  of  the  turbulent  mean  state  squared  Brunt-Vaisala  (buoyancy) 
frequency  (N2).  A  dramatic  decrease  of  N2  in  the  jet  core,  through  turbulent  dynamics  and 
mixing,  is  clearly  noticed  (it  takes  negative  values  in  case  3  for  some  vertical  levels  implying 
the  existence  of  unstably  stratified  layers).  A  N2  profile  of  this  kind,  with  a  localized  mini¬ 
mum  in  the  mixing  region,  is  sometimes  referred  to  as  possessing  a  lN2  notch’.  It  is  known 
to  be  a  mean  state  configuration  favoring  emission  of  gravity  waves  from  tropospheric  jet 
streams  ([Lott,  Kelder  k  Teitelbaum  (1992)];  [Sutherland  k  Peltier  (1995)]).  We  note  that 
although  the  basic  state  Brunt-Vaisala  profile  has  only  a  smooth  doubling  at  the  jet  lev¬ 
els,  the  time-varying  mean  state,  in  our  DNS,  evolves  to  a  profile  having  such  a  N 2  notch. 
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[Nastrom  &  Eaton  (1997)]  found  a  localized  decrease  of  N2  at  the  tropopause  level  in  several 
winter  seasonal  profiles,  which  adds  credence  to  our  simulation  results  (cf.  their  Fig.  5b). 
The  vertical  profiles  of  the  gradient  Richardson  number  Rig  =  N2q/  (dU/dzf  at  the  quasi¬ 
equilibrium  are  shown  in  figure  10b  for  cases  1,  2  and  3.  Here  N2q  =  Nd{d)  jdz  +  N2K2(z) 
are  the  normalized  Brunt-Vaisala  frequency  profiles  at  quasi-equilibrium.  The  gradient 
Richardson  number  Rig  profile  at  quasi-equilibrium  is  very  similar  to  those  obtained  by 
[Bedard,  Canavero  &  Einaudi  (1986)]  as  shown  in  figure  3a  and  to  the  profiles  for  cases  4,  5 
and  6  in  figure  3b.  Comparing  cases  1,  2  and  3,  we  can  see  that  an  increase  of  stratification 
causes  the  curves  of  Rig  to  move  closer  to  the  core.  On  the  other  hand,  increase  or  decrease 
in  the  forcing  term  also  changes  the  vertical  position  of  the  Rig  curves,  which  is  clear  when 
we  compare  cases  1  and  5. 

Visualizations  of  the  flow  fields  are  shown  in  figures  11-14  for  a  better  understanding 
of  the  flow  structure.  The  velocity  vectors,  fluctuation  temperature,  spanwise  vorticity  and 
local  gradient  Richardson  number  on  a  typical  vertical  plane  cross-section  (where  the  hori¬ 
zontal  scale  is  taken  to  be  ax)  are  plotted  here  for  cases  1  and  2.  As  can  be  seen  in  figure  11a, 
Kelvin-Helmholtz  type  instability  at  the  edges  of  the  jet  produces  counterclockwise  flow  on 
top  and  clockwise  flow  on  the  bottom.  This  agrees  with  the  contour  levels  of  the  spanwise 
vorticity  in  figure  12a.  Entraining  fluid  from  above  and  below  the  jet  generates  fingers  of 
hot  and  cold  fluid,  as  evident  from  temperature  contours  shown  in  figure  lib.  The  temper¬ 
ature  fluctuation  at  the  center  is  comparatively  small.  Figure  12b  shows  the  local  gradient 
Richardson  number  at  the  same  vertical  plane  cross-section.  High  values  are  obtained  away 
from  the  core  while  values  close  to  the  center  are  small.  The  action  of  gravity  waves  prop¬ 
agating  away  from  the  jet  stream  leads  to  enhanced  thermal  fluctuations  at  other  vertical 
levels.  Peaks  of  temperature  variance  (see  also  figure  18)  are  shifted  relative  to  peaks  of 
velocity  variances  shown  in  figures  15-16.  These  phenomena  have  been  observed  in  the  at¬ 
mospheric  measurements  by  [Nastrom,  Gage  &  Ecklund  (1986)]. 

The  fields  for  case  2  are  plotted  in  figures  13  and  14  (stronger  background  stratification 
compare  to  case  1).  The  vortices  and  the  waviness  are  observed  in  the  vector  projections 
(figure  13a).  The  corresponding  spanwise  vorticity  field  shown  in  figure  14a  has  high  inten¬ 
sity  zones  concentrated  around  the  center.  The  most  significant  differences  between  cases  1 
and  2  are  the  amplitude  of  temperature  fluctuations,  which  are  intensified  in  case  2.  The 
maxima  and  minima  of  temperature  fluctuations  in  case  2  are  closer  to  the  core;  that  is 
za«2  instead  of  za  «  5  in  case  1  (cf.  figure  18a).  The  local  gradient  Richardson  number 
shown  in  figure  14b  is  not  significantly  different  from  case  1. 
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The  variances  of  the  streamwise  velocity,  (u'2)/U(0)2  are  shown  in  figure  15.  Cases  1,  3 
and  5  have  similar  magnitudes  while  case  2  has  the  smallest  normalized  variance.  In  general, 
the  variances  show  two  peaks  located  roughly  at  the  edges  of  the  jet,  at  za  sa  ±0.70  where 
dUb/dz  is  largest.  This  agrees  with  observations  of  [Beland  (1993)]  that  the  turbulence  in¬ 
tensity  is  highest  at  the  edges  of  the  jet  stream.  The  value  of  (u2)/U{0)2  is  the  lowest  in 
case  2,  and  decays  more  rapidly  away  from  the  core  due  to  the  higher  stratification.  In 
contrast,  the  values  for  cases  1  and  3  are  similar,  suggesting  that  stratification  in  these  cases 
is  too  low  to  have  any  significant  effect  on  velocity  variances.  This  is  also  true  for  cases 
4  &  5.  The  variances  of  the  spanwise  velocity  (v'2)/U( 0)2  are  shown  in  figure  16.  Their 
vertical  profiles  exhibit  two  peaks,  but  their  magnitudes  are  generally  much  smaller  than  for 
the  corresponding  streamwise  components.  For  cases  4  &  5,  the  values  of  spanwise  velocity 
variance  at  the  core  increase  as  the  Reynolds  number  increases.  Typically,  the  spanwise 
component  has  smaller  relative  value  in  the  TKE  budget  than  the  streamwise  component, 
but,  as  the  Reynolds  number  increases  intercomponent  energy  transfer  also  increases,  thus 
reducing  anisotropy  of  the  field.  Vertical  velocity  variances  (w'2)/U(0)2  are  presented  in 
figure  17 .  They  show  a  single  peak  at  the  center  of  the  jet  core.  Again,  case  2  is  significantly 
smaller  in  value  than  cases  1  and  3,  in  view  of  increased  stratification.  The  magnitudes  for 
cases  4  and  5  are  roughly  the  same  as  in  cases  1  and  3. 

The  normalized  temperature  variances  (#2)/Z7(0)2  shown  in  figure  18  indicate  that  tem¬ 
perature  and  vertical  velocity  tend  to  be  out  of  phase,  the  former  peaks  where  the  latter 
drops  to  small  values.  The  asymmetry  in  all  cases  are  very  prominent,  a  result  of  the  jump 
in  Brunt- Vaisala  frequency  profile  at  the  center.  As  stratification  decreases,  the  turbulence 
penetrates  deeper  into  the  surrounding  fluid,  spreading  the  peak  of  temperature  variance 
further.  In  case  2,  the  upper  peak  of  (i9'2)/[7(0)2  is  located  at  za  «  2.4;  in  case  1,  at  za  »  5 
and  in  case  3  at  za  8  which  are  even  further  away  from  the  core  than  the  peak  of  <  -d  >. 
The  vertical  levels  of  peaks  for  <t9'2)/f7(0)2  tend  to  vary  inversely  with  the  strength  of  strat¬ 
ification.  For  cases  4  &  5,  the  temperature  variance  peaks  roughly  at  the  same  vertical 
levels  as  in  case  1.  However,  the  magnitude  of  the  peaks  decreases  as  the  Reynolds  number 
increases. 

The  normalized  horizontal  heat  fluxes,  (u'#)/U(0)2,  are  plotted  in  figure  19.  They  are 
negative  immediately  above  the  center  and  positive  below  the  center.  This  is  most  prominent 
in  case  2.  This  is  in  agreement  with  the  previously  discussed  figures  11  and  13  where  fingers 
of  hot  fluid  are  entrained  into  the  core  from  above  and  fingers  of  cold  fluid  from  below.  Fur¬ 
ther  away  from  the  jet,  Id'u')  changes  sign  again.  This  is  caused  by  the  counterflow  which 
carries  cooler  fluid  to  the  top  and  hotter  fluid  to  the  bottom,  resulting  in  a  positive  and  a 
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negative  secondary  peak  above  and  below,  respectively. 


The  vertical  heat  fluxes  (w'lf) /U(Q)2  are  plotted  in  figure  20.  For  cases  1,  2  and  3,  they 
are  all  negative,  confirming  that  the  heat  is  moving  from  hotter  (lighter)  upper  domain  to 
cooler  (heavier)  lower  domain  (down-gradient).  Counter-gradient  heat  flux  ((it/#)  positive) 
has  been  observed  in  others  studies  only  when  the  gradient  Richardson  number  is  relatively 
large  (Rig  «  0.5  -  1.0,  [Gerz,  Schumann  k  Elghobashi  (1989)]).  For  case  3,  the  heat  flux 
is  almost  zero  at  the  center  corresponding  to  the  near  neutral  stratification.  The  equation 
(2.28)  for  mean  temperature  suggests  that  the  mean  temperature  is  related  to  the  gradient 
of  the  vertical  heat  flux.  This  agrees  with  the  fact  that  the  peak  of  (i9)  in  each  case  shown  in 
figure  9  corresponds  to  the  maximum  rate  of  change  of  (it/#).  As  a  result,  the  peak  of  (i9) 
is  situated  further  away  than  the  peak  of  («/'#).  In  figure  20,  the  profiles  of  (it/#)/F(0)2 
for  cases  1,  2  and  3  are  significantly  different.  If  the  profiles  are  not  normalized  (not  shown 
here),  case  2  has  the  largest  peak  at  the  core  of  the  flow.  When  stratification  decreases,  the 
amount  of  heat  flux  at  the  core  of  the  flow  decreases,  until  it  drops  nearly  to  zero  for  case  3. 
At  the  same  time,  heat  flux  at  the  edges  increases  and  generates  two  peaks  at  za  &  7  and 
—6.5  in  case  3.  The  doubling  of  the  Brunt-Vaisala  profile  results  in  the  asymmetry  of  the 
heat  flux  profiles. 

Anisotropy  in  stably  stratified  mixing  layers  has  been  studied  by  [Smyth  &  Mourn  (2000)]. 
They  found  that  the  vorticity  becomes  isotropic  as  the  Reynolds  number  increases.  Veloc¬ 
ity  has  a  similar  trend  but  it  is  possible  that  in  real  atmospheric  situations  some  degree 
of  anisotropy  will  be  retained.  Indeed,  anisotropy  of  turbulent  fields  observed  in  strato¬ 
spheric  measurements  is  discussed  in  [Vinnichenko  (1980)].  Recent  observations  of  turbu¬ 
lence  around  a  winter  jet  stream  using  the  ‘Egrett’  meteorological  research  aircraft  actually 
show  that  anisotropy  persists  at  the  smallest  scale  observed  (~  Ira) 

([Cote,  Roadcap,  Wroblewski,  Dobosy  k  Crawford  (2003)]). 

4  Variability  of  length  scales  and  gradient  Richardson 
number 

The  behaviour  of  stably  stratified  turbulent  flows  can  be  characterized  in  terms  of  length 
scales  where  buoyancy,  shear,  inertial  and  viscous  influences  become  dominant  or  length- 
scale  ratios  that  indicate  the  relative  magnitudes  of  these  influences.  Some  commonly  used 
length-scales  are  the  shear  length  scale  Ls,  Ellison  scale  Le,  buoyancy  scale  Lb  and  Ozmidov 
scale  L0  defined  in  table  3.  Here  e  is  the  dissipation  rate  of  turbulent  kinetic  energy  intro- 
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duced  in  eq.  (2.26).  Relationships  between  these  length  scales  have  been  investigated  exten¬ 
sively,  for  example,  in  the  atmospheric  boundary  layer  by  [Hunt,  Kaimal  &  Gaynor  (1985)], 
through  atmospheric  radar  observations  by  [Eaton  &  Nastrom  (1998)],  in  numerical  sim¬ 
ulations  by  [Itsweire,  Koseff,  Briggs  &  Ferziger  (1993)])  and  using  laboratory  experiments 
([Rohr,  Itsweire,  Helland  &  Van  Atta  (1988)];  [De  Silva  &  Fernando  (1992)]  (see  [Fernando  (2002)] 
for  recent  discussions).  From  these  scales,  several  important  scale  ratios  can  be  constructed. 

For  example,  the  ratio  Rig  =  L\jL\  already  shown  in  figure  3b  and  figure  10b  is  the  gradient 
Richardson  number  Rig. 


Table  3:  Definition  of  length  scales. 


streamwise 

spanwise 

vertical 

other 

Kolmogorov  scale,  L2 

- 

- 

- 

Shear  scale,  L2 

{u2) /(dU/dz)2 

- 

(w2)/ (dU/dz)2 

q2/{dUjdz)2 

Ellison  scale,  L2 

- 

! 

- 

0 »'2>/  (§Y 

Buoyancy  scale,  L\ 

- 

- 

(®'2 )/K 

i2/K 

Ozmidov  scale,  L2 

- 

- 

Body  forces,  such  as  background  buoyancy,  typically  tend  to  influence  motion  scales 
that  are  larger  than  certain  ‘outer’  scales  (Ozmidov,  buoyancy,  Ellison  and  shear).  The 
scales  of  active  Kolmogorov  turbulence  are  those  scales  which  are  greater  than  the  Kol¬ 
mogorov  scale  but  smaller  than  any  outer  scale.  [Ozmidov  (1965)]  suggested  that  buoy¬ 
ancy  effects  are  important  on  a  scale  proportional  to  L0  defined  in  table  3.  The  Ozmi¬ 
dov  scale  is  regarded  as  the  outer  scale  above  which  the  fluid  turbulence  is  significantly 
affected  by  buoyancy  (cf.  [Phillips  (1972)],  (1991)).  Also,  the  dynamics  of  a  stably  strat¬ 
ified  turbulent  flow  depend  on  the  rate  at  which  fluid  elements  mix  with  each  other  and 
change  their  density  in  the  interior  of  the  flow.  [Ellison  (1957)]  introduced  the  scale  Le 
now  known  as  the  overturning  or  Ellison  scale  which  uses  the  density  (temperature)  fluctua¬ 
tions  as  an  indicator  of  the  vertical  scale  turbulent  motions;  Le  is  a  typical  vertical  distance 
traveled  by  fluid  particles  before  either  returning  towards  their  equilibrium  level  or  mixing 
([Rohr,  Itsweire,  Helland  &  Van  Atta  (1988)]), 

The  ratio  Ls/L b  =  RilJ2  accounts  for  the  relative  scales  of  influence  of  shear  and  buoy¬ 
ancy.  If  Rig  >  1,  then  the  buoyancy  influence  occurs  at  a  smaller  scale  than  the  shear,  and  if 
active  turbulence  exists,  then  e  ~  (w'2)^2/Lb  and  hence  Lb  ~  (w'2)^2/N  ~  (eLh)l^/N  ~  L0 
indicating  that  the  Ozmidov  scale  La  represents  the  scale  beyond  which  the  turbulent  eddies 
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are  influenced  by  stratification  ([Fernando  &  Hunt  (1996)]).  Note  that,  in  general,  L0  should 
be  contrasted  from  the  buoyancy  scale  Lb  at  which  the  fluid  parcel  displacements  (either  due 
to  turbulence  or  waves)  are  constrained  by  the  buoyancy  forces.  The  Kolmogorov  scale  Lv 
characterizes  the  smallest  scale  of  active  turbulence  and  L0  ~  Ln  then  characterizes  the  com¬ 
plete  suppression  of  active  turbulence  by  stratification,  leading  to  states  with  significant  wave 
activity.  The  Ellison  scale  Le  is  representative  of  the  vertical  motions  of  a  stratified  turbulent 
fluid,  in  that  it  specifies  how  fluctuations  associated  with  such  motions  (due  to  eddies  or  waves 
or  combination  thereof)  are  related  to  the  vertical  scale  of  motions.  Given  that  both  Lb  and 
Le  respond  to  the  presence  of  waves  and  turbulence  alike,  [Hunt,  Kaimal  &  Gaynor  (1985)] 
have  suggested  that  the  ratio  Le/Lb,  which  is  called  the  temperature  fluctuation  parameter, 
is  insensitive  to  the  presence  of  waves,  with  an  estimated  value  of  Le/Lb  =  1.0  ±  0.5,  For 
the  weakly  stable  case  of  Rig  <  1  with  intense  turbulence,  the  Ozmidov  scale  has  a  limited 
significance  and  large  L0JLb  values  are  possible. 

In  figures  21a,b  three  length  scales  L0,  Le  and  Lv  are  compared  for  cases  1  and  2,  re¬ 
spectively.  The  scales  are  normalized  by  1/a,  and  for  the  sake  of  clarity  their  logarithmic 
values  are  shown;  note  the  difference  in  vertical  variability  of  different  length  scales.  While 
Le  and  L0  are  bell  shaped  and  peak  roughly  at  the  jet  core,  Lv  becomes  minimum  therein. 
Ozmidov  and  Kolmogorov  scales  differ  significantly  in  the  core  indicating  active  turbulence 
with  strong  mixing,  concurrent  with  the  low  Rig  values  in  the  core  (figure  10b).  The  ratio  of 
L0/ Ln  decreases  markedly  along  the  jet  edges,  indicating  increased  inhibition  of  mixing  and 
impact  of  buoyancy  on  much  smaller  scales.  Beyond  za  «  ±6,  the  turbulent  fluctuations 
die  down  (figures  15-18),  so  as  the  vertical  heat  flux  (figure  20).  The  scales  Le  and  L0  in  the 
upper  domain  drop  faster  than  in  the  lower  domain  because  of  higher  stratification.  Figure 
21b  is  a  plot  of  length  scales  for  case  2,  which  is  generally  similar  to  case  1.  The  Kolmogorov 
scale  remains  unchanged,  but  two  peaks  of  L0  can  be  seen  in  the  jet  shear  layers.  However, 
the  levels  where  Lv  «  L0  and  Ln  ~  Le  now  occurs  at  a  much  lower  vertical  levels  za  &  ± 3, 
a  signature  of  decreased  spatial  influence  of  jet  due  to  increased  stratification.  These  ob¬ 
servations  are  consistent  with  the  fields  shown  in  figures  11-14.  The  decades  of  separation 
between  Ln  and  the  outer  scales  in  the  zone  of  active  turbulence  do  guarantee  the  effective 
resolution  of  all  scales.  For  ease  2,  figure  21b,  Le  becomes  greater  than  L0  for  za  >  2  (the 
ratio  Le/L0  crosses  unity  at  za  «  2);  from  figure  10b  these  are  also  the  vertical  levels  where 
Rig  >  0.25,  with  pronounced  temperature  fluctuations  and  very  patchy  turbulence  in  that 
region.  This  is  in  general  agreement  with  figure  10  of  [Schumann  &  Gerz  (1995)]  who  have 
plotted  ratios  of  Lt/L0  vs  Richardson  number  and  observed  Le/L0  <  1  for  Rig  small,  and 
Le/L0  1  at  Rig  ~  0,25  (in  the  context  of  uniform  linear  vertical  shear  and  stratification). 
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Figures  22  a,b  show  the  variations  of  L0,  Lb  and  Ls  for  cases  1  and  2,  respectively.  In 
figure  22a,  within  the  jet  core  \za\  <  2  (except  at  the  center)  Ls  «  Lb,  indicating  a  low 
Rig(«  1)  conditions  prone  for  instabilities  and  turbulence.  In  the  core  La  »  Ls,  and 
somewhat  larger  than  Lb  indicating  that  turbulence  in  this  region  is  active  and  intense. 
Note  that  Ls  <  Lb  implies  that  the  vertical  fluid  particle  displacements  are  sufficiently  large 
to  be  deformed  by  the  mean  shear,  promoting  the  interaction  between  Reynolds  stresses  and 
mean  shear  to  yield  sustained  turbulence  in  the  core  (\za\  <  2),  which  is  clear  from  figures 
15-19.  In  figure  22b  which  has  higher  stratification  than  case  1  we  have  L0  <  Lb  throughout 
the  vertical  profile,  indicating  prominent  wave  activity. 

Since  Rig  is  small  (<  0,25)  in  the  jet  core  (except  the  center  za  =  0),  the  turbulence 
therein  is  intense  and  the  temperature  is  well  mixed,  leading  to  low  temperature  fluctu¬ 
ation  levels  in  the  core  (figures  3b,  lib,  18).  In  the  higher  Rig  region  where  buoyancy- 
dominated  turbulence  prevails,  temperature  fluctuations  become  maximum,  thus  providing 
conditions  for  higher  optical  turbulence.  This  is  a  region  of  interest  where  the  length  scale 
ratio  L0/Lb  <  1  points  to  significant  wave  activity  (figure  24a).  Detailed  energy  budgets 
(see  Section  5;  figures  25-26)  show  that  in  this  region  the  shear  production  is  vanishing  and 
the  energy  is  received  via  energy  transfer  from  the  jet  core  region.  However,  the  energy 
deposited  in  this  region  by  the  pressure-velocity  term  is  much  smaller  (and  sometime  neg¬ 
ative)  than  that  gained  by  the  nonlinear  terms,  the  wave  activity  appears  to  be  dominated 
by  non-linear  waves.  Given  that  the  wave  energy  does  not  propagate  beyond  \za\  «  6,  it 
appears  that  these  non-linear  waves  break  and  cause  intermittent  turbulent  patches  in  the 
region  (2  <  \za\  <  6),  which  is  consistent  with  the  high  local  Rig  regions  of  figures  12b  and 
14b.  Such  breaking  events  can  sustain  some  heat  flux  (as  shown  in  figures  19-20)  in  this 
layer  while  maintaining  high  temperature  fluctuations. 

At  the  edges  of  the  jet  (case  1,  figure  22a),  in  a  thin  layer  around  \za\  «  1.8,  we  have 
Ls  >  La,Lb  indicating  that  the  turbulence  production  by  shear  and  Reynolds  stress  inter¬ 
action  is  impeded.  Rig  in  this  region  is  large  (figure  3b)  and  hence  the  active  turbulence  is 
possible  when  Rig  locally  drops  below  a  critical  value  (0.25)  to  produce  Kelvin-Helmholtz 
type  instabilities.  Figures  11-14  show  that  indeed  the  turbulence  is  intermittent  at  the  jet 
edges  and  Kelvin-Helmholtz  type  instabilities  are  present.  At  larger  distances  from  the 
jet  center  2.2  <  \za\  <  4,  Lb  ~  Ls,  indicating  marginal  conditions  for  the  maintenance  of 
stratified  turbulence;  figure  12  shows  that  turbulence  in  this  region  is  very  patchy.  Also,  pro¬ 
nounced  temperature  fluctuations  in  this  region  (with  Le  >  L0  for  case  2)  points  to  possible 
wave  activity.  For  \za\  >  4,  the  scales  L0  and  Lb  are  markedly  smaller  than  Ls,  confirming 
the  lack  of  sustained  turbulence  when  shear  is  neglegible.  This  is  generally  consistent  with 
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figures  11-14  and  velocity  variances  shown  in  figures  15-17. 


Other  length  scale  ratios  are  also  of  interest.  As  stated  above  [Hunt,  Kaimal  &  Gaynor  (1985)] 
predicted  Le/Lb  ratio  of  unity  in  stable  atmospheric  boundary  layers,  and  found  a  ratio  of 
about  0.8  in  their  field  measurements.  The  ratios  of  Le/Lb  in  our  DNS  for  cases  1,  2,  4  and 
5  are  plotted  in  figure  23a.  The  ratio  is  smallest  around  the  center  of  the  jet  and  increase 
at  the  edges.  It  drops  to  0.1  at  the  center  for  case  1.  In  the  core,  Le/Lb  for  case  2  (with 
strong  stratification)  is  higher  («  0.25)  than  that  of  case  1  indicating  that  turbulence  is 
more  buoyancy  influenced.  At  the  center,  the  ratios  for  cases  4  and  5  are  even  smaller  than 
case  1.  For  all  cases,  the  ratios  increase  close  to  unity  away  from  the  core,  indicating  the  in¬ 
creasing  importance  of  wave-turbulence  activities  at  such  vertical  levels.  In  case  2,  the  ratio 
approaches  a  peak  value  of  unity  at  around  za  ±4,  while  for  cases  1,  4  and  5,  the  peaks  of 
Le/Lb  are  observed  at  za  «  ±6.  From  figure  18,  these  peaks  approximately  correspond  to 
the  observed  peaks  of  (i 9'2 )  for  both  cases  1  and  2.  Higher  stratification  pushes  the  vertical 
peak  levels  of  the  ratio  closer  to  the  core. 

There  are  other  methods  to  determine  the  generation  of  waves  from  jet  turbulence.  For 
example,  using  linear  equations,  [Stewart  (1969)]  suggested  that  the  phase  spectra  of  w  is  90° 
ahead  of  d  if  the  motion  is  wave-like.  The  method  was  used  by  [Lienhard  &  Van  Atta  (1990)] 
but  they  did  not  obtain  conclusive  results.  Another  method  is  the  ‘penetration  condition’ 
derived  and  used  by  [Sutherland  &  Peltier  (1993)],  (1995)  also  from  linearized  equations.  In 
the  two-dimensional  context,  they  used  the  method  to  demonstrate  the  generation  of  internal 
gravity  waves  from  a  jet  profile  very  similar  to  ours.  The  Brunt- Vaisala  profile  was  twice 
larger  in  the  upper  domain  than  in  the  lower  domain  with  a  notch  at  the  center.  In  their 
simulations,  the  absolute  magnitude  of  penetration  ratio  is  greater  than  unity,  the  peak  of 
the  Reynolds  stress  moves  up  the  domain  as  time  evolves,  suggesting  the  presence  of  linear 
gravity  waves.  If  the  penetration  ratio  were  smaller  than  unity,  the  peak  should  also  travel 
upward,  but  they  find  a  significantly  smaller  magnitude.  They  attribute  this  discrepancy  to 
nonlinear  effects.  The  simulations  were  done  in  a  2D  decaying  context.  In  our  simulations, 
we  cannot  reproduce  the  upward  motion  of  the  Reynolds  stress  peaks  (in  quasi-equilibrium 
state)  and  it  is  not  certain  whether  similar  results  can  be  obtained  in  a  3D  context.  Never¬ 
theless,  our  results,  which  use  vertical  variability  of  LejLb  as  an  indicator  of  wave  activity, 
corroborate  their  findings.  Remarkably,  the  curves  for  Le/Lb  show  little  variation  for  cases 
4,5  and  6  (figure  23b),  and  saturate  with  increased  resolution;  vertical  resolution  doubles 
from  cases  4  and  5  (512  vertical  levels)  to  case  6  (1024  vertical  levels).  This  guarantees  that 
our  resolution  of  the  outer  scales  is  grid  independent  and  physically  significant. 
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Figure  24a  shows  the  ratio  of  L0/Lb,  which  is  an  indicator  of  the  levels  of  turbulence  and 
waves.  When  L0jLb  <  1,  the  waves  activity  is  considered  prominent  over  turbulence  and 
vice  versa.  Note  how  this  ratio  changes  over  the  vertical  spread  of  the  jet.  The  value  of 
L0/Lb  ~  1  observed  outside  the  jet  core  in  the  shear  layer  is  of  the  same  order  of  magnitude 
as  that  measured  by  [Hebert,  Mourn,  Paulson  k  Caldwell  (1992)]  in  a  turbulent  patch  in 
Pacific  equatorial  (jet-like)  undercurrent. 

The  ratio  ( L0/Le )2(/3  is  sometimes  called  the  turbulent  Froude  number,  Frt  and  is  an 
alternative  to  the  bulk  Richardson  number  ([Shih,  Koseff,  Ferziger  k  Rehmann  (2000)]). 
Frt  for  cases  1,  2,  4  and  5  are  plotted  in  figure  24b.  The  curves  are  bell  shaped,  with 
cases  4  and  5  having  the  largest  values,  followed  by  case  1  and  then  case  2.  The  val¬ 
ues  at  the  core  are  all  larger  than  unity,  but  Frt  drops  to  less  than  unity  in  regions 
where  the  stratification  is  dominant.  In  laboratory  grid-generated,  salt-stratified  turbulence, 
[Rohr,  Itsweire,  Helland  k  Van  Atta  (1988)]  have  measured  the  ratio  Le/L0  from  which  the 
corresponding  turbulent  Froude  numbers  Frt  can  be  deduced.  For  N  =  0.96,  they  obtain 
values  from  0.1  to  1.8  as  the  downstream  position  increases,  which  correspond  to  ( L0/Le )2/3 
values  of  4.6  to  0.68.  Their  values  are  in  the  same  range  as  those  depicted  in  figure  24b. 


5  Budgets 

Budget  terms  were  computed  in  the  quasi-equilibrium  state  in  order  to  characterize  turbu¬ 
lent  transport  processes  occurring  at  various  vertical  levels.  Our  budget  equations  presented 
in  this  section  are  based  on  the  decomposition  of  velocity,  temperature  and  pressure  fields 
into  the  mean  (horizontally  averaged)  and  fluctuating  components  u',  v' ,  w',  9'  and  p' .  This 
decomposition  is  given  by  eqs.  (2.6)-(2.10).  The  budget  equations  (5.29)-(5.32)  are  obtained 
after  substitution  of  the  decomposition  (2.6)-(2.10)  into  the  governing  equations  and  per¬ 
forming  operations  of  horizontal  averaging  ( ).  More  discussions  of  the  budget  equations 
based  on  the  decomposition  (2.6)-(2.10)  can  be  found  in  [Stull  (1988)].  We  report  the  bud¬ 
get  terms  of  (u2),  (v'2),  (w'2)  and  <t9'2)  for  case  2  (high  stratification),  where  the  heat  flux 
is  larger  compared  to  other  terms. 

The  budget  equation  for  (u2)  is  : 

d 

—{u2)=  +  11^1 2)  +  IUtf  2)  +  IVtf  2)  (5.29) 

where 


V2)  =  —2  v(Vu'-Vu'), 
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TT  ,du' . 

"<*'*>  =  2<p&>' 

"V»>  =  |(^((x'2»  -  (^»  , 

'V>  =  -2(u'w')^U. 

The  budget  terms  are  normalized  by  the  quantity  Z7(0)3a;.  The  normalized  terms  for  (u2) 
are  plotted  in  figure  25a  for  the  strongly  stratified  case  2.  The  only  source  term  is  the  shear 
production  (IV{u>2)).  It  has  two  peaks,  located  on  either  side  of  the  model  jet  stream,  in 
agreement  with  the  graph  of  { u  2)  in  figure  15.  The  energy  produced  is  either  dissipated  or 
transferred  to  other  components  through  the  dissipation  (//u'2^)  and  pressure  redistribution 
terms.  The  pressure  term  also  consists  of  two  peaks  directly  opposite  to  the  peaks 
of  the  production  terms.  This  term,  is  effective  in  transferring  energy  from  the  streamwise 
direction  and  depositing  it  to  the  spanwise  and  vertical  directions.  The  dissipation  term  has 
the  largest  values  at  the  center,  and  it  decreases  away  from  the  center.  The  transport  term 
moves  energy  from  the  peaks  to  the  center  and  a  lesser  amount  to  the  outer  edges 

of  the  jet. 

The  budget  equation  for  the  spanwise  velocity  is  : 
d 

—(v2)=  I{v>2)  +  IIW2)  +  +  IVtf 2) ,  (5.30) 

where 


I{v' 2>  =  -2is(W  •  W) , 

=  ^(^((»2))  -  (w'tx'2))) , 

IV^  =  -2<»V)|V. 

The  budget  terms  for  ( v  2)  are  shown  in  figure  25b.  The  source  is  the  pressure  (i7(i/2)) 
which  transfers  energy  from  the  streamwise  to  the  spanwise  direction.  The  magnitude  of  the 
pressure  term  is  small  compared  to  the  similar  term  in  the  budget  equation  for  the  vertical 
velocity  (figure  26a).  The  dissipation  term  (/(l/2)),  in  contrast  to  the  streamwise  component, 
has  two  side  peaks.  The  transport  term  (. Ill ^)  moves  energy  from  either  side  to  the  center 
and  away  from  the  jet.  The  shear  production  term  IV^v>2^  in  the  spanwise  direction  is  small 
compared  to  other  terms. 
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(5.31) 


The  budget  equation  for  the  vertical  velocity  is: 
d 

0^(W  2)  =  I(w'2)  +  +  ln{w'2)  +  IV(w'2)  > 

where 

1{w'2)  =  -2  v{Vw'-Vw'), 

uv*)  =  -  <»'3>  -  2<pv: », 

=  Sw. 

The  budget  terms  for  (u/2)  are  shown  in  figure  26a.  The  pressure  redistribution  term 
(^(w'2))  takes  a  larger  portion  of  energy  from  the  streamwise  direction  than  the  spanwise 
term  II(v>2y  The  energy  is  then  either  dissipated  (i)w'2p  or  converted  to  potential  energy 
through  the  buoyancy  term  {IV{w>2)).  The  dissipation  term  here  has  magnitude  similar  to 
that  in  the  spanwise  direction;  but  this  term  is  now  dominant  at  the  center.  The  transport 
term  (III{w'2))  plays  a  more  significant  role  than  in  the  spanwise  direction,  in  that  it  trans¬ 
ports  energy  from  the  jet  shear  layer  area  to  the  outer  region  (which  may  be  responsible 
for  the  intermittent  wave  breaking  activity  seen  in  figures  13-14).  The  buoyancy  term  is 
negative  and  heat  flux  is  thus  down  gradient  (from  hot  to  cold).  No  counter  gradient  heat 
flux  was  observed  from  the  equilibrium  statistics  (figure  20). 

The  budget  equation  for  the  temperature  is  : 

d 

2)  =  J^2)  4-  II^>  2}  +  1 1  I^i  2^ ,  (5.32) 

where 

1{0'2)  =  -2/c(W  •  W'> , 

=  -2(wW)  (NK2  +  = -2(w'l>')^ , 

"V->  = 

The  budget  terms  for  (d"1)  are  shown  in  figure  26b.  The  gradient  production  term 
y)  is  the  only  source  in  the  balance  of  potential  energy.  It  has  a  peak  above  za  fa  2 
where  the  buoyancy  flux  for  case  3  is  larger.  The  dissipation  term  (J^)  has  roughly 
the  same  shape  as  the  production  term  but,  of  course,  with  the  opposite  sign.  This  may 
explain  why  countergradient  fluxes  are  not  observed  in  figure  26a;  temperature  anomalies 
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dissipate  in  the  generation  region,  without  producing  restratification  effects  that  would  lead 
to  countergradient  fluxes.  The  transport  term  (III^i2\)  is  relatively  small. 

The  normalized  shear  production  (IV^)  in  the  streamwise  direction  for  cases  1,  2  and 
3  is  shown  in  figure  27a.  Cases  1  and  3  have  nearly  identical  shear  production.  Case  2, 
which  is  more  stably  stratified,  has  smaller  shear  production  at  each  vertical  level  compare 
to  cases  1  and  3.  Thus,  the  normalized  shear  production  decreases  with  increasing  stratifica¬ 
tion  in  our  flows.  This  is  consistent  with  the  conclusion  in  [Holt,  Koseff  k  Ferziger  (1992)] 
that  stratification  does  not  directly  reduce  the  growth  of  q 2  (i.e.  not  a  kinetic  energy 
sink).  Rather,  the  effect  of  stratification  is  indirect  in  suppressing  the  production  (see 
also  [Rohr,  Itsweire,  Helland  k  Van  Atta  (1988)]).  The  normalized  total  dissipation  (J^/2)  + 
I(v> 2)  +I{W'2))  for  the  three  cases  is  shown  in  figure  27b.  The  dissipation  for  case  2  decreases 
following  the  decrease  in  production.  The  ratio  of  production  over  dissipation  is  shown  in  fig¬ 
ure  27c.  The  three  curves  nearly  collapse  into  one,  suggesting  that  production  and  dissipation 
maintain  the  same  ratio  regardless  of  stratification.  However,  the  ratio  changes  with  za  and 
peaks  at  the  edges  of  the  jet,  where  shearing  is  maximum,  za  ml.  Around  the  jet  production 
and  dissipation  are  not  at  equilibrium  and  the  transport  must  be  taken  into  account.  Another 
relevant  parameter  is  the  growth  parameter  F  as  defined  in  [Holt,  Koseff  k  Ferziger  (1992)] 


production  —  buoyancy  ! 

l/V^I 

1  - 1 

IV(w'2) 

dissipation 

'2> 

1  +  IV2) 

l  +  l 

(5.33) 


The  dependence  of  F  on  za  is  shown  in  figure  27d.  It  varies  roughly  between  0  and  2.5 
depending  on  the  vertical  level.  We  note  that  F  is  roughly  equal  to  unity  in  homogeneous 
simulations  for  critical  Richardson  number,  e.g.  [Holt,  Koseff  k  Ferziger  (1992)],  and  in 
experiments,  e.g.  [Strang  k  Fernando  (2001)],  suggesting  that  TKE  transport  plays  a  more 
significant  role  in  the  situations  considered  herein. 


6  Conclusions 


A  spectral  domain  decomposition  method  was  developed  to  study  inhomogeneous  shear- 
stratified  turbulent  flows  in  the  ABL  context  with  focus  on  the  detailed  structure  of  a  turbu¬ 
lent  jet  in  the  tropopause  region.  In  contrast  to  most  previous  work,  the  method  employed 
here  permits  more  realistic  boundary  conditions  in  the  vertical  direction  and  allows  shear 
and  stratification  profiles  to  adjust  during  flow  evolution.  Quasi-equilibrium  solutions  were 
obtained  via  long-time  integration  of  governing  equations.  In  this  state,  the  micro-scale 
Reynolds  number  and  other  averaged  quantities  fluctuate  within  5  %  of  their  mean  value. 
In  previous  numerical  studies,  the  micro-scale  Reynolds  number  was  temporally  evolving 
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except  at  a  certain  “stationary”  Richardson  number. 

The  quasi-equilibrium  jet  produced  during  the  simulations  consists  of  a  jet  core  sur¬ 
rounded  by  shear  layers,  beyond  which  lay  a  (return)  counterflow.  The  Richardson  number 
within  the  jet  core  is  low,  the  Ozmidov  scale  is  high  and  turbulent  energy  production  peaks 
within  the  core,  thus  producing  sustained  turbulence.  The  Richardson  number  increases 
beyond  the  jet  shear  layers  in  such  a  way  that  the  shear  length-scale  exceeds  the  buoyancy 
scale,  thus  impeding  sustained  energy  transfer  to  turbulence  by  the  mean  flow.  The  turbu¬ 
lence  in  this  region  is  patchy.  Beyond  the  jet  edges  are  regions  with  moderate  Richardson 
number  characterized  by  small  Ozmidov  to  Buoyancy  scale  ratio,  indicating  the  presence  of 
stably  stratified  turbulence;  the  shear  production  rate  therein  is  vanishingly  small  and  the 
energy  supply  to  this  region  occurs  via  propagating  non-linear  gravity  waves  (the  non-linear 
transport  term  in  the  energy  budget  exceeds  the  pressure  transport  term).  This  is  supported 
by  a  clear  difference  between  the  core  and  the  edges  of  the  jet  with  respect  to  the  ratios  of 
length  scales.  For  example,  the  ratio  Le/Lh  is  significantly  smaller  than  unity  at  the  cen¬ 
ter  while  approaching  unity  for  all  cases  away  from  the  jet  core,  implying  the  increasing 
importance  of  wave  activity  there.  This  indicates  that  edge  regions  above  and  below  the 
jet  stream  maximum  shear  are  the  regions  where  much  of  the  strong  interactions  between 
the  mean  flow,  locally  generated  Kelvin-Helmholtz  instabilities,  turbulence  and  propagating 
gravity  wave  instabilities  take  place.  This  is  important  for  ABL  optical  propagation  codes.  It 
is  also  in  agreement  with  the  observational  studies  in  [Bedard,  Canavero  k  Einaudi  (1986)] 
and  [Dalaudier,  Sidi,  Crochet  k  Vernin  (1994)]. 

Theoretical  analysis  of  stratified  shear  flows  has  further  established  that  propagating 
wave  modes  are  forced  at  the  shear  layers  by  a  Kelvin-Helmholtz  type  shear  instability;  this 
shear  instability  then  supports  wave  radiation  when  the  density  stratification  outside  of  the 
shear  layer  is  sufficiently  large  ([Lott,  Kelder  k  Teitelbaum  (1992)]);  this  is  borne  out  by  the 
present  simulations.  Our  DNS  further  show  that  with  decreasing  background  stratification, 
the  vertical  levels  of  enhanced  nonlinear  wave  interactions  significantly  separate  away  from 
the  levels  of  peak  shear  production.  Such  vertical  variability  of  turbulence  statistics  must 
be  factored  into  ABL  propagation  codes  and  parameterizations  of  C. % 

Our  simulations  were  performed  with  lower  Reynolds  numbers  than  in  the  atmosphere 
and  some  of  our  results  are  effected  by  the  low  Reynolds  number  simulations.  However,  cer¬ 
tain  important  properties  of  computed  turbulent  fields  are  not  only  in  qualitative  but  also 
in  quantitative  agreement  with  observed  geophysical  flows.  For  example,  computed  ratios 
and  variability  of  turbulence  natural  length  scales  (Ozmidov,  buoyancy,  shear  and  Ellison)  are 
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broadly  consistent  with  the  observational  data  presented  in  [Bedard,  Canavero  &  Einaudi  (1986)], 
[Eaton  &  Nastrom  (1998)],  [Nastrom,  Gage  &  Ecklund  (1986)],  [Hunt,  Kaimal  &  Gaynor  (1985)] 
and  [Hebert,  Mourn,  Paulson  &  Caldwell  (1992)].  The  ratios  of  turbulence  outer  scales  are 
found  to  be  in  agreement  with  observations  and  they  do  not  change  with  increasing  Reynolds 
number  and  saturate  with  numerical  resolution.  This  is  demonstrated  in  figure  3b  for  the 
ratio  of  shear  and  buoyancy  length  scales  where  vertical  profile  of  Rig  =  L\jL\  is  shown  for 
cases  4,  5  and  6  in  figure  3b;  and  the  ratios  Le/Lj,  for  cases  4,  5  and  6  in  figure  23b. 

Description  of  turbulence  and  especially  turbulence  closure  models  in  stratospheric  layers 
in  the  ABL  context  require  parameterizations  and  scaling  laws  for  the  length  scales.  As 
there  are  a  number  of  different  natural  length  scales  for  stably  stratified  turbulent  flows  with 
shear,  it  is  not  clear  without  additional  theoretical  arguments,  numerical  or  experimental 
data  which  scales  are  appropriate  for  the  present  case.  One  of  the  main  objectives  of  this 
research  is  to  establish  for  non-homogeneous  stratified  shear  flows  in  the  ABL  context  in 
the  vicinity  of  the  tropopause  the  length  scales  over  which  different  phenomena  occur,  their 
relative  orders  of  magnitude  and  variability.  In  future  work  we  will  investigate  in  more  depth 
the  parameterizations  for  various  turbulence  statistics  in  the  context  of  nonhomogeneous 
tropopausal  turbulence  in  the  ABL  context  to  support  the  ABL-ADA  (atmospheric  decision 
aid). 
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Figure  1:  (a)  An  example  of  a  jet  stream  profile  from  Bedard,  Canavero  & 
Einaudi  (1986)  (reproduced  with  permission  from  the  American  Meteorological 
Society).  The  dotted  line  is  obtained  from  observation  while  the  solid  line  is  a 
numerical  fit.  (b)  The  normalized  velocity  profiles  E4  at  the  basic  state  of  the 
present  computations.  The  vertical  axis  is  za  where  a  -  16,  and  -8.0  <  za  < 
8.0.  Only  the  middle  portion  of  the  computational  box  is  shown  here. 


Figure  2:  Brunt-Vaisala  frequency  profile  recorded  in  the  atmosphere  by  Beland 
(1996)  with  balloon  measurements;  vertical  axis  -  height  (km)  and  horizontal 
axis  -  Brunt-Vaisala  frequency  (1/s). 
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Figure  3:  (a)  An  example  of  a  gradient  Richardson  number  profile  correspond- 
ing  to  figure  1(a),  measured  by  Bedard,  Canavero  &  Einaudi  (1986)(repro- 
dueed  with  permission  from  the  American  Meteorological  Society);  (b)  Gradient 
Richardson  number  profiles  at  quasi-equilibrium  from  the  numerical  simulations 
for  cases  4,  5  (2562  x  512  resolution)  and  case  6  (5122  x  1024  resolution). 
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Figure  4:  The  Reynolds  number  based  on  the  length  scale  Ld  and  q  for  5122  x 
1024  resolution,  case  6,  Red  =  qLd/v. 
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Figure  5:  (a)  Distribution  of  collocation  points  in  primary  and  auxiliary  subdo¬ 
mains.  The  vertical  lines  separate  the  subdomains,  (b)  Distribution  of  primary 
subdomain  boundaries  in  the  vertical  direction. 


Figure  7:  The  time  evolution  of  vertically  integrated  values  of  (u2)  (I),  (v2) 
(II),  {w  2)  (III),  (i?  2)  (IV)  and  | d{u  2)/dt\  (V);  {*)  denotes  horizontal  averaging. 
Horizontal  axis  is  given  in  non-dimensional  time  Sbt. 


Figure  8:  Normalized  mean  (horizontally  averaged)  velocity  in  the  streamwise 
direction  for  cases  listed  in  Table  1,  horizontal  axis  -  U(z)/U(0)  where  U  = 
Ub  +  {«),  vertical  axis  -  za  :  (a)  cases  1  (solid),  2  (dash)  and  3  (dash-dot)  ;  (b) 
cases  4  (solid)  ,  5  (dash)  and  6  (dash-dot). 


Figure  9:  Normalized  mean  (horizontally  averaged)  temperature,  ($)(z)/ 17(0), 
in  the  quasi-equilibrium  state  for  cases  listed  in  Table  1,  horizontal  axis  -  (i9), 
vertical  axis  -  za;  (a)  cases  1  (solid),  2  (dash)  and  3  (dash-dot)  ;  (b)  cases  4 
(solid),  5  (dash)  and  6  (dash-dot). 


Figure  10:  (a)  The  square  of  the  Brunt-Vaisala  frequency  profile,  N2K2,  and  (b) 
Gradient  Richardson  number  profiles  at  quasi-equilibrium  from  the  numerical 
simulations  for  cases  1  (solid),  2  (dash)  and  3  (dash-dot). 


Figure  13:  (a)  Instantaneous  velocity  vectors  and  (b)  fluctuating  temperature 
component  on  a  vertical  plane  for  case  2.  The  horizontal  axis  is  xa. 


Figure  14:  (a)  Instantaneous  spanwise  vorticity  and  (b)  local  gradient  Richard¬ 
son  number  on  a  vertical  plane  for  case  2.  The  horizontal  axis  is  xa. 
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Figure  15:  The  normalized  variances  of  the  streamwise  velocity,  (u2)  / 77(0) 2  as  a 
function  of  the  normalized  height :  (a)  cases  1  (solid),  2  (dash)  and  3  (dash-dot); 
(b)  cases  4  (solid)  and  5  (dash). 


Figure  16:  The  normalized  variances  of  the  spanwise  velocity,  (v'2)/U(0)2  :  (a) 
cases  1  (solid),  2  (dash)  and  3  (dash-dot)  ;  (b)  cases  4  (solid)  and  5  (dash). 
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Figure  17:  The  normalized  variances  of  the  vertical  velocity,  {w'2)/U( 0)2  :  (a 
cases  1  (solid),  2  (dash)  and  3  (dash-dot)  ;  (b)  cases  4  (solid)  and  5  (dash). 


Figure  18:  The  normalized  variances  of  the  temperature,  {d'2)/U(Q)2  :  (a)  case 
1  (solid),  2  (dash)  and  3  (dash-dot)  ;  (b)  cases  4  (solid)  and  5  (dash). 


Figure  19:  The  normalized  profile  of  the  horizontal  heat  flux,  (■d'u')/U(0)2  :  (a 
cases  1  (solid),  2  (dash)  and  3  (dash-dot)  ;  (b)  cases  4  (solid)  and  5  (dash). 


Figure  20:  The  normalized  profile  of  the  vertical  heat  flux,  (#w')/?7(0)2  :  (a 
cases  1  (solid),  2  (dash)  and  3  (dash-dot)  ;  (b)  cases  4  (solid)  and  5  (dash). 
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Figure  25:  (a)  The  normalized  horizontally  averaged  budgets  in  the  equation 
for  the  streamwise  velocity  variance  for  case  2.  Labels  in  the  figure  correspond 
to  I  -  dissipation,  II  -  pressure,  III  -  transport,  and  IV  -  shear  production, 
(b)  The  normalized  budgets  in  the  equation  for  the  spanwise  velocity  variance. 
Labels  in  the  figure  correspond  to  /  -  dissipation,  II  -  pressure,  III  -  transport 
and  IV  -  shear  production  (term  IV  is  nearly  zero). 


Figure  26:  (a)  The  normalized  horizontally  averaged  budgets  in  the  equation 
for  the  vertical  velocity  variance  for  case  2.  Labels  in  the  figure  correspond  to 
I  -  dissipation,  II  -  pressure,  III  -  transport,  and  IV  -  buoyancy,  (b)  The 
normalized  budgets  in  the  equation  for  temperature  variance.  Labels  in  the 
figure  correspond  to  I  -  dissipation,  II  -  production  and  III  -  transport. 


:malized  dissipation;  (c) 
meter  for  cases  1  (solid), 


e  nearly  identical. 


